Comprehensive analysis of LILR family genes expression and tumour‐infiltrating immune cells in early‐stage pancreatic ductal adenocarcinoma

Abstract Leucocyte immunoglobulin‐like receptors (LILRs) are closely related to tumourigenesis, but their clinical value in early‐stage pancreatic ductal adenocarcinoma (PDAC) after pancreaticoduodenectomy remains unknown. Kaplan–Meier and Cox proportional hazards regression models is used to investigate the association between LILR expression and prognosis in tumour biopsies and peripheral blood mononuclear cells. Risk score was calculated for each patient based on the prognostic model. DAVID, STRING, GeneMANIA, and GSEA were used to conduct pathway and functional analyses. The CIBERSORT algorithm is used to analyse tumour‐infiltrating immune cells. Survival analysis showed that high levels of LILRA4 (p = 0.006) and LILRB4 (p = 0.04) were significantly associated with better overall survival. High levels of LILRA2 (p = 0.008) and LILRB4 (p = 0.038) were significantly associated with better relapse‐free survival. JAK‐STAT signalling pathway, regulation of T cell activation, regulation of the immune effector process, and tumour necrosis factor superfamily cytokine production were involved in molecular mechanisms that affected poor prognoses in the high‐risk group in GSEA. CIBERSORT demonstrated that the high‐risk group had significantly higher infiltrating fraction of memory‐activated CD4 T cells and activated NK cells and lower fraction of resting dendritic cells and neutrophils. LILRB4 plays crucial roles in affecting the clinical outcomes of early‐stage PDAC.

pancreatic cancer has risen from 2.5% in 1970 to 9% in 2019, there is still a large gap in the improvement of the clinical outcomes of other tumour types, such as breast and colorectal cancers [3][4][5]. Lack of biological indicators for early effective screening strategies and novel therapies in PDAC are the main reasons for this discrepancy [6,7]. Current studies have shown that tumour markers, such as various commonly used carcinoembryonic antigens and carbohydrate antigen 199, are not effective in diagnosing early PDAC due to their low sensitivity [7]. However, early detection strategies involving endoscopic ultrasound and magnetic resonance imaging have not been validated in randomised trials in high-risk PDAC patients [4]. Therefore, identifying effective screening indicators and treatment-related targets is crucial to improve the prognosis of pancreatic cancer.
PDAC is similar to most adenocarcinomas, with a massive fibrotic stroma, which plays an important role in the local inflammatory microenvironment of tumours [8][9][10]. The microenvironment of PDAC consists of abundant deposition of the extracellular matrix, low vessel density, cancerassociated fibroblasts, and immune/inflammatory cells, which are closely related to tumour growth and progression and the infiltration of immune cells, similar to other tumours [2,[11][12][13][14][15]. Immune cells may be found in solid tumours, and different types of immune cells have different effects on the clinical prognosis of tumours [16,17]. In particular, increased number of dendritic cells (DCs) is associated with improved prognosis in various types of human cancers. DC maturation is a prognostic indicator; moreover, chronic inflammation and presence of M2 macrophages facilitate tumour growth and spread [18,19].
Leucocyte immunoglobulin-like receptors (LILRs) are a family of receptors with extracellular immunoglobulin domains; the gene coding region is located on the chromosome region 19q13.4, also known as CD85, ILT, and LIR, which have immunomodulatory effects on a variety of immune cells [20,21]. LILRs include subfamily A and subfamily B (LILRA1-6 and LIRB1-5, respectively). The LILR receptor in the A subfamily is an activation receptor that contains tyrosinebased immune receptor activation motifs (ITAMs), while the LIR receptor in the B subfamily contains multiple tyrosinebased cytoplasmic immune receptor inhibition motifs (ITIMs) [22,23]. LILRs are related to the human killer cell inhibitory receptor (KIR) family. Both have similar Ig-like structures and cytoplasmic signal domains. Although the expression of KIR is limited to natural killer (NK) cells, LILRs are found in a variety of cells, including NK, T lymphocytes, B lymphocytes, and myeloid cells (monocytes, macrophages, dendritic cells, and granulocytes [24]). The transmembrane domain of the LILRA receptor contains a charged arginine or lysine residue associated with the FcR ligand containing (ITAXI/Lx 6-12 YxxI/L) ITAM [25]). ITAM activation recruits Syk/ZAP70 family kinases to drive downstream activation pathways, which are important for immunity [26]. By contrast, the LILRB receptor contains the cytoplasmic (S/I/ V/LxYxxI/V/L) ITIM domain, which recruits the phosphatase SHP1/SHP2/SHIP containing the Src homology 2 domains, thereby inhibiting the immune signalling cascade. SHP/SHIP phosphatase activity is essential to maintain immune homoeostasis [27]. LILRAs and LILRBs belong to the LILRs group, and the main function of LILRAs involves immune activation, while the role of LILRBs is immune suppression [28][29][30]. A study on oestrogen receptor-positive breast cancer reported that after receiving neoadjuvant endocrine therapy, compared with patients with low expression of the LILRA2 gene in the tumour biopsies, patients with high LILRA2 expression showed significant tumour shrinkage, which was beneficial to breastconserving surgery [31]. Lu et al. found that the Semaphorin-4A gene stimulated the CD4+ T cells and regulated Th2 T cell differentiation by binding to LILRA2, thereby initiating an immune response [32]. A previous study reported that LILRA4 was also expressed in cancer cells, which was associated with impairment of plasmacytoid dendritic cells (pDCs) in the microenvironment of cancers [33]. A genomewide analysis of copy number variations in ovarian cancer has shown that duplicate mutations in LILRA6 were associated with susceptibility to high-grade serous ovarian cancer [34].
Human Leucocyte Antigen-G (HLA-G) positive expression in gastric cancer patients indicated a poor prognosis, and its possible mechanism may be that HLA-G combined with LILRB1 inhibited NK cell proliferation and function [35]. A series of studies have shown that in a variety of tumours, such as colon cancer, non-small cell lung cancer, and hepatocellular carcinoma, overexpressed LILRB2 was associated with a poor prognosis [36][37][38][39][40][41]. LILRB4, which is also an immunosuppressive receptor, is similar to LILRB3, so it can interact with HLA-G to inhibit the activation of immune cells, which mainly include NK and T cells [42,43]. High LILRB4 expression has also been reported to be associated with tumour progression and poor prognosis. LILRB4 was associated with impaired T cell responses in pancreatic cancer, and antagonistic LILRB4 was the key to successful immunotherapy [44,45]. In leukaemia, tumour cells disable immune checkpoint blockade therapy through the LILRB4 signalling, and blocking LILRB4 can prevent the development and metastasis of tumour cells. The potential mechanism is that LILRB4 changes the tumour microenvironment, resulting in immune suppression [46,47]. In colorectal cancer, the LILRB4 gene was found to be highly expressed in cancer tissue, and the expression level of the LILRB4 gene was negatively correlated with the density of CD45RO + T cells in the cancer tissue, and high LILRB4 gene expression was a biomarker of poor clinical prognosis for colorectal cancer [48]. Moreover, in gastric cancer, LILRB4 was significantly related to the pathological grade [49]. Therefore, we know that LILR family genes play an important role in tumourigenesis, development, and clinical prognosis.
However, the effect of LILRs on clinical outcomes in patients with PDAC remains unknown. This study aimed to investigate the association between LILRs gene expression and prognosis and reveal possible mechanisms in pathway enrichment and tumour immune cell infiltration in early-stage PDAC.

| Data mining and processing
The gene expression profiles and clinical information of PDAC were obtained from The Cancer Genome Atlas (TCGA, data release 21.0, December 10, 2019) and normalised by the 'DESeq' and 'Limma' package in R (version 3.6.1; www-project.org) [50]. The inclusion criteria for cases in this study were as follows: (i) patients with pancreaticoduodenectomy and pathologically confirmed as PDAC, (ii) according to the Seventh American Joint Committee on Cancer, patients with postoperative specimen pathological stage were stage I or II, and pathology stages I and II were defined as early-stage PDAC, and (iii) patients with complete clinical prognosis data. To further characterise the expression levels of LILRs between PDAC patients and healthy control patients in the peripheral blood mononuclear cells (PBMCs), GSE74629 and GSE49641 were obtained from the gene expression omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo) database. GSE55643 was used to verify differences in gene expression levels. Data from the GEO database were normalised by the 'normalizeBetweenArrays' function of 'Limma' package in R.

| Bioinformatics and correlation analysis of LILR genes
First, we analysed the differences in LILR gene expression between pancreatic cancer and adjacent tissues and constructed violin plots using the 'wilcox.test' function of 'vioplot' package in R. We then used the 'corrplot' package in R to conduct Pearson's correlation analyses between LILR genes and then constructed corresponding correlation plots. We used the online tool DAVID Bioinformatics Resources (version 6.8, https://david.ncifcrf.gov) to conduct functional annotation of gene ontology (GO) terms and enrichment of the Kyoto Encyclopaedia of Genes and Genomes (KEGG) for all LILR genes. All 10 genes from the LILR gene family were included in this analysis. [51,52]. Finally, we analysed the interrelationships of gene-gene and protein-protein interactions between LILRs genes using online analysis tools: GeneMA-NIA( [53]) (http://genemania.org) and STRING (https:// string-db.org) [54].

| Survival analysis and prognostic model construction
Using the median expression of each LILR gene as the cut-off point, we defined the high-and low-expression groups. Kaplan-Meier analysis with the log-rank test was used as a univariate analysis to assess the association between clinical factors, gene expression, and clinical outcomes, including OS and relapse-free survival (RFS). Clinical factors with p-values of the log-rank test less than 0.05 were used as corrective factors when constructing the Cox proportional hazards regression model of each LILR gene. Then, the genes that were significantly associated with clinical outcomes were recruited to construct a prognostic model. The formula of risk scores was as follows: where β is the regression coefficient from the multivariate Cox proportional regression model of individual genes. We used median risk scores as the cut-off point and divided the patients with PDAC into higher and lower risk groups. A receiver operating characteristic (ROC) curve was used to evaluate the predictive power of the prognostic model, which was conducted using the 'roc' function of 'ROC' package in R [55]. Subsequently, we combined the risk score with clinical factors of p-values less than 0.05 for a nomogram analysis, which was also completed using the 'rms' and 'survival' packages of R [56]. Afterwards, we analysed differences in expressions of prognosis-related genes between PDAC patients and healthy control patients in the peripheral blood and PBMCs and visualised the results using the 'ggplot' function and t-test of 'ggplot2' package in R. ROC curves were constructed to show the diagnostic efficiency of prognosis-related genes in the diagnosis of PDAC by using the 'pROC' package in R.
To verify the stability of the model, we used the k-fold method that was used for cross-validation of the model, 'caret' package to create cross-validation data set, and 'survi-valROC' package to calculate AUC values.

| GSEA
After constructing the prognostic model, we found that compared with the lower risk group, the higher risk group had a worse clinical prognosis. To investigate the potential molecular mechanism of adverse prognosis in the higher risk group, we performed GSEA between both groups [57]. The Molecular Signatures database (MSigDB) C2 and C5 gene sets were used in the GSEA. The C2 gene set mainly included KEGG and related signalling pathway analysis, while the C5 gene set mainly included GO enrichment analysis. The selection criteria for statistically significant gene sets were a false discovery rate (FDR) less than 0.25 and a p-value less than 0.05.

| Analysis of tumour-infiltrating immune cells
LILR genes were closely related to the regulation of immune cells, so we performed a tumour-infiltrating immune cell (TIIC) analysis between the higher risk group and lower risk group. The CIBERSORT algorithm was used to perform tumour-infiltrating immune cell analysis, which is a gene expression-based deconvolution algorithm, which can measure 22 types of characteristic immune cell compositions in RNA GAO ET AL. mixtures from many tissues, including solid tumours [58]. The algorithm based on the mRNA expression sequence was conducted using R, and the set parameter was 1000 permutations.

| Statistical analyses
All statistical analyses and plots were performed using SPSS 23.0 and R (version 3.6.1). A p-value of less than 0.05 was considered statistically significant. Association analyses between clinical factors and gene expression and clinical outcomes were tested using Kaplan-Meier analysis and a Cox proportional hazard regression model, and hazard ratios (HRs) with 95% confidence intervals were used to describe the relative risks. A two-sided t-test was used to evaluate the differential expression analysis of LILR gene expression and the percentage of tumour-infiltrating cells between the two groups. In GSEA, the Benjamini-Hochberg method was used to adjust multiple tests, which defined the meaning of FDR; in our study, FDR <0.25 was considered statistically significant.

| Data processing
After standardised processing, 112 patients with PDAC and 31,777 genes who met the inclusion criteria from TCGA were enroled in our study. Clinical factors included age, sex, history of chronic pancreatitis and alcohol, tumour size, pathological stage, neoplasm histological grade, targeted molecular therapy, radiation therapy, residual resection, and information of clinical outcomes (Table S1). Overall survival data were available for all samples in the 112 PDAC patients, but 19 cases were missing in the relapse-free survival data.

| Bioinformatics and correlation analysis of LILR genes
We obtained 10 LILR gene expressions from the pancreatic cancer expression profile of TCGA. These included LILRA1, LILRA2, LILRA4, LILRA5, LILRA6, LILRB1, LILRB2, LILRB3, LILRB4, and LILRB5 (Table S3). We performed a differential expression analysis of those genes between pancreatic cancer tissues and adjacent tissues. The results showed that compared with adjacent tissues, LILRA1, LILRA2, LILRA4, LILRA6, LILRB1, LILRB2, LILRB3, and LILRB4 had a statistically significant lower expression in cancer tissues (p < 0.05; Figure 1a). We then performed correlation analyses on these 10 LILR genes and visualised the results with correlation diagrams as shown in Figure 1b. The results showed that they were all positive correlations, and most of them were strong correlations (correlation coefficient ≥0.6). Some of them were very strong correlations (correlation coefficient ≥0.8), such as LILRB1 and LILRB4, LILRB2, and LILRB3.
Furthermore, we conducted functional annotation and pathway enrichment for the LILR genes. As shown in the results of Figure 2a and Table S2, the LILR genes were mainly involved in the adaptive immune response, the regulation of immune response, the immune system process, leucocyte differentiation, MHC class I protein binding, transmembrane signalling receptor activity, and inhibitory MHC class I receptor activity. Notably, the Fc receptor-mediated inhibitory signalling pathway, immune response-regulating cell surface receptor signalling pathway, regulation of immune response, regulation of immune system process, and immune response were in the same regulatory network, which suggested that the LILRs genes played an important role in the immune response and regulation. The results of protein-protein interaction (PPI) analysis showed that LILR genes were closely associated with each other, and they were related to the HLA family and integrin family genes ( Figure 2b). The results from GeneMA-NIA were similar, showing a close relationship between LILR genes ( Figure 2c).

| Survival analysis and prognostic model construction
In the analysis of the association between the above clinical factors and clinical outcomes, we found that four variables (neoplasm histological grade, targeted molecular therapy, radiation therapy, and residual resection) were significantly associated with OS, and two factors (neoplasm histological grade and residual resection) were significantly associated with RFS as shown in Table S1. Clinical factors of appeal were included in the Cox proportional hazard regression model as corrective factors. Other factors were not significantly associated with clinical outcomes.
We then conducted a joint analysis of genes that were significantly associated with clinical outcomes based on gene expressions. According to the level of LILRA4 and LILRB4  Table 3 and Figure 4b).
Subsequently, we constructed a prognostic model for OS based on LILRA4 and LILRB4, and a prognostic model for RFS was based on LILRA2 and LILRB4. In this study, the β value of the Cox proportional hazard regression model was negative; so to facilitate analysis and graph display, we performed a logarithmic transformation of raw risk scores. The specific formula of the prognostic model was as follows: risk scores (OS) = ( −log10(LILRA4 expression � −0.673 + LILRB4 � −0.26)))+ 4, risk scores (RFS) = (-log10-(LILRA2 expression � −0.123 + LILRB4 � −0.095)])) + 4. Here, the constant 4 was used to make the final risk scores a positive number. Based on median risk scores as the cut-off point, patients were divided into higher risk and lower risk groups. From the results of Table 4 and Figure 5a,b, we found that the higher risk score group was significantly associated with poor OS compared to the lower risk score group (p = 0.039, HR = 1.75, 95% CI: 1.03-3.00). In the prognostic model of RFS, the result showed that patients with PDAC after pancreaticoduodenectomy with higher risk scores suggested a poor RFS (p = 0.007, HR = 2.90, 95% CI: 1.33-6.15) (Table 5; Figure 6a

46
-Analysis of the tumour biopsies revealed that LILRA2, LILRA4, and LILRB4 were associated with the prognosis of PDAC patients. To test if these genes had prognostic potential in PBMCs as well, GSE49641 and GSE74629 datasets were analysed. We found that the difference in the expression of these genes in cancer patients and the healthy controls was not statistically significant in GSE49641 ( Figure S2A-C Supporting Information), but LILRA2 and LILRB4 were significantly higher in the peripheral blood of cancer patients in GSE74629 ( Figure S2D-F Supporting Information). Due to significant differences in the expression of LILRA2 and LILRB4 in GSE74629 in peripheral blood of PDAC patients, we decided to construct an AUC curve to evaluate the diagnostic efficacy of LILRA2 and LILRB4 in the diagnosis of PDAC. The AUCs of LILRA2 and LILRB4 were 0.853 and 0.896, respectively ( Figure S2G in Supporting Information), which showed the strong ability of LILRA2 and LILRB4 in diagnosing PDAC by detecting the expression levels in peripheral blood. Then, we set to confirm whether the lower expression of these genes is in the tumour tissue compared to the normal peri-tumoral pancreatic tissues holds in a validation dataset, namely GSE55643. LILRA2, LILRA4, and LILRB4 showed low expressions in tumour tissues, in which LILRB4 was significantly lowly expressed (p < 0.001; Figure S3A in Supporting Information).

| GSEA
GSEA was used to determine the potential molecular mechanisms that affected poor prognoses in the higher risk score groups. Our results indicated that has04650 (natural killer cellmediated cytotoxicity), hsa04062 (chemokine signalling), and hsa04630 (JAK-STAT signalling pathway) were enriched in the C2 gene set, which may be involved in affecting OS and RFS (Figure 7a). In the C5 gene set, GO 0050863 (regulation of T cell activation), GO 0002697 (regulation of immune effector process), and GO 0071706 (tumour necrosis factor superfamily cytokine production) were enriched and significantly correlated with OS and RFS (Figure 7b).

| Analysis of TIICs
According to the previous analysis from the TCGA database, LILRs may influence postoperative clinical outcomes of patients with PDAC by participating in the regulation of the immune response and the function of immune cells, but whether the immune infiltration contributed to it was unknown. Hence, we used CIBERSORT to evaluate the differential percentages of TIICs between the higher and lower risk score groups. CIBERSORT is a deconvolution algorithm based on gene expression; it is combined with 22 leucocyte gene signature matrices, which are a defined 'barcode' with 547 gene expression signatures to distinguish the subgroups of 22 type immune cells. Figure 8a,b shows the landscape of tumourinfiltrating immune cells for 112 early-stage PDAC cases. Tumour-infiltrating immune cell populations for each immune cell type are shown in Figure 8c. We then performed correlation analyses on the relative percentages of these 22 immune cells. The result showed that there was a strong correlation between naïve CD4 T cells and memory B cells (correlation coefficient = 0.73) (Figure 9). Furthermore, the results of differential percentages of the tumour-infiltrating immune cells showed that there was a difference in T cell CD8 (p = 0.02), memory-activated CD4 T cells (p = 0.001), activated NK cells (p = 0.011), M1 macrophages (p = 0.025), resting dendritic cells (p = 0.024), eosinophils (p = 0.02), and neutrophils (p = 0.018) between both groups using the OS prognostic model (Figure 10a). In addition, the percentages of tumourinfiltrating cells in the RFS prognostic model, including plasma cells (p = 0.038), memory-activated CD4 T cells (p = 0.008), gamma delta T cells (p = 0.023), activated NK cells (p = 0.008), M2 macrophages (p = 0.039), resting dendritic cells (p = 0.004), and neutrophils (p = 0.002), were statistically significant between the two groups ( Figure 10b). In addition, in order to understand the TIICs in PMBCs and solid tumour tissues, we analysed GSE74629 and found that the results were inconsistent with the analysis results in TCGA, possibly because the cell types of the components were different ( Figure S3D in Supporting Information).

| DISCUSSION
In this study, we investigated the association between LILR genes and clinical prognosis and diagnosis in early-stage PDAC. Our results suggested that compared with adjacent tissues, LILRA1, LILRA2, LILRA4, LILRA6, LILRB1, LILRB2, LILRB3, and LILRB4 were significantly overexpressed in cancer tissues. Multivariate analyses of the Cox proportional hazard regression model showed that higher expressions of LILRA4 and LILRB4 were significantly associated with better OS, and there was a significant association between lower expression of LILRA2 and LILRB4 and a worse RFS. However, inconsistent with the results of other tumours, high expressions of LILRB2 and LILRB4 were associated with better clinical outcomes in our study, while other studies showed that patients with high expressions of T A B L E 2 Combined survival analysis of LILRA4 and LILRB4 gene expression with overall survival in early-stage PDAC.  LILRB2 and LILRB4 predicted worse clinical outcomes [45][46][47][48][49][50][51][52]. The inconsistency of the results was confusing and interesting for further study. Tumours are involved in an abnormally complex process, and the entire process cannot be explained by individual gene events. Therefore, we conducted a joint analysis and constructed a prognostic model based on gene expressions significantly associated with clinical outcomes. The joint effect of clinical variables and LILR expression indicated that patients with two risk factors had higher hazard ratios than those with only one risk factor. Our findings T A B L E 4 The association analysis between risk scores and overall survival in early-stage PDAC. suggested that the risk score can be used to evaluate the clinical outcomes of patients with early-stage PDAC. The AUC on the ROC curves was slightly smaller than some well-known prognostic scores, such as the Glasgow prognostic score, or the modified Glasgow prognostic score. However, the HR value was similar [59,60]. The difference was that the prognostic score of this study was also effective in the evaluation of RFS. Moreover, based on the prognostic model, tumourinfiltrating immune cells were investigated. Our findings can therefore help assess the RFS for PDAC patients and identify immunotherapeutic targets. The prognostic model included LILRA and LILRB genes, consistent with real-world studies, and were grouped according to the level of risk score to investigate the possible mechanism of risk score influencing clinical prognosis using GSEA analysis, which may make the results more reliable. GSEA results revealed that the following pathways may be involved in regulating the potential molecular mechanisms that affect clinical prognosis; the JAK-STAT signalling pathway, regulation of T cell activation, regulation of immune effector processes, and tumour necrosis factor superfamily cytokine production were enriched and significantly correlated with the OS and RFS. Numerous studies have shown that the activation of the JAK-STAT signalling pathway promotes the development and progression of tumours, including pancreatic cancer [61][62][63][64]. Therefore, we have a reason to speculate that in the low-risk group, activation of the JAK-STAT signalling pathway played an important role in affecting the OS and RFS. At the same time, many studies have reported that the JAK-STAT signalling pathway is closely related to immune evasion, immune regulation, immune cell differentiation, and drug resistance [64][65][66].

Groups
The results of GSEA also showed enrichment of the immune-related signalling pathway, which suggested that the JAK-STAT signalling pathway seems to have some association. Hence, we investigated the differential percentages of TIICs between the higher and lower risk score groups. Our results found that in a prognostic model of OS and RFS, there were significantly different percentages in memory-activated CD4 T cells, activated NK cells, resting dendritic cells, and neutrophils between the two groups. A series of studies have shown that TIICs are closely related to tumour development, drug response, and clinical prognoses in a variety of tumours [67][68][69][70][71][72]. In these studies, a high fraction of infiltrating and activated CD4 F I G U R E 7 The results of GSEA for the OS (a) and RFS (b) prognostic model. The median expression of each LILR gene was used as the cut-off point, which defined the higher and lower risk groups for the early-stage PDAC. GAO ET AL. memory T cells, resting dendritic cells, and neutrophils indicated a good clinical prognosis, which is consistent with previous studies. However, in contrast to previous studies, the high-risk group had high percentages of activated NK cells, while the low-risk group had lower fractions of activated NK cells, which was similar to the relationship between the expression of LILRB4 and clinical prognosis in early-stage PDAC. In general, NK cells play an anti-tumour role and are regarded as prognostic factors for a good prognosis in some tumours, but tumourinfiltrating NK cells may also play a role in promoting tumours [73]. NK cells are divided into two subtypes according to the cell surface antigens CD56 and CD16; namely CD56dim/ CD16+ and CD56bright/CD16-. CD56dim/CD16+ NK cells have a high cytotoxic potential, while the CD56bright/ CD16-subtype has a low cytotoxic potential but can secrete cytokines that promote angiogenesis. The vast majority of tumour-infiltrating NK cells are the CD56bright/CD16-subtype and have been reported to have impaired cytotoxicity and cytokine-producing functions and promote angiogenesis, thus playing a role in promoting tumours [74,75]. Therefore, it is possible that in the high-risk group, a high fraction of tumourinfiltrating NK cells plays a role in promoting tumours, which is a reasonable explanation. The infiltrating immune cells are part of the tumour microenvironment. In the prognosis model, the combined effect of the infiltrating immune cells indicates the possible mechanism of poor prognosis in the high-risk group. The bioinformatics analysis in solid tissues was based on a complex statistical method and was calculated using a formula.
CIBERSORT was verified in lung specimens obtained during surgical resection of early-stage non-small-cell lung carcinomas and disaggregated lymph node biopsies from follicular lymphoma by flow cytometry. The results of CIBERSORT were significantly correlated with flow cytometry measurements ( [76]), but the real correlation in PDAC needs to be verified by additional studies.

| CONCLUSION
In this study, we investigated the associations between LILRs genes and clinical prognosis in early-stage PDAC. The results revealed that LILRA4 and LILRB4 were significantly F I G U R E 9 Correlogram of Pearson's correlation analysis of tumour-infiltrating immune cells; red denotes a positive correlation, blue denotes a negative correlation, and the shade of colour represents the size of the correlation coefficient. associated with OS, and LILRA2 and LILRB4 were associated with the RFS. LILRB4 was significantly related to the prognosis of PDAC patients, but the prognostic values of LILRA2 and LILRA4 need further validation. The prognostic model suggested that patients with early-stage PDAC with higher risk scores had worse clinical outcomes. The results of potential molecular mechanistic analyses suggested that in the prognostic model, the JAK-STAT signalling pathway, immunerelated signalling pathways, and tumour-infiltrating immune cells may play a crucial role in affecting the clinical outcomes, but further functional experiments are needed for confirmation. F I G U R E 1 0 A violin diagram of differential infiltrating infraction analysis for the OS prognostic model (a) and RFS prognostic model (b).